This document contains the analyses of the manuscript entitled “Regional heterothermy is not associated with active heat dissipation by the horns in the Rhinoceros Beetle (Megasoma gyas)” by Danilo Giacometti, Luiz Henrique Lima Silva, Guilherme Gomes, José Eduardo de Carvalho, and Alexandre V. Palaoro.
To maintain thermal balance, animals typically rely on the control of physiological and behavioral processes. Some animals, however, possess exaggerated morphological structures that aid in thermoregulation, particularly in the dissipation of excess heat when body temperatures rise. Although widespread in the animal kingdom, the role of animal weapons—exaggerated morphological structures with multiple characteristics—has received little attention in thermoregulation studies. Here, we investigated whether the horns of the Rhinoceros Beetle (Megasoma gyas) acted as a thermal window. We heated live and dead beetles to 30 ºC and allowed them to cool to 20 ºC, while measuring surface temperature changes in four body regions: the cephalic horn, the thoracic horn, the elytra, and the abdomen. If the horn actively dissipated heat, the cephalic horn would show the lowest cooling rate among body regions. Contrary to this expectation, we found that the cephalic horn had the highest cooling rate, followed by the abdomen, thoracic horn, and elytra, respectively. This suggests that the horns are not used actively for heat dissipation in M. gyas. Live specimens had lower cooling rates than control specimens across all body regions, which could be an indication of hemolymph flow modulating heat dissipation. However, whether beetles show circulatory responses to temperature is currently unknown. The low cooling rate of the elytra can be explained by the high density of flight muscles in the thorax. Apart from their role in heat generation, flight muscles may aid in heat dissipation by pumping hemolymph across tagmata or through the low-insulated cuticle to prevent thoracic overheating. Our study also demonstrates that invertebrates show regional heterothermy even in instances unrelated to exercise or stress. We propose that regional heterothermy results from both active (control of hemolymph flow) and passive (heat dissipation through poorly insulated structures) processes within individual.
ggtheme <- function(base_size=12, base_line=0.3) {
theme(
text = element_text(size=base_size),
line = element_line(size=base_line, linetype="solid"),
axis.text.x = element_text(size=base_size*0.8, colour='black', hjust=0.5, vjust=1, angle=0),
axis.text.y = element_text(size=base_size*0.8, colour='black', hjust=1, vjust=0.5, angle=0),
axis.line.x = element_line(size=base_line),
axis.line.y = element_line(size=base_line),
axis.title.x = element_text(size = base_size, vjust = 1, margin=unit(c(3,0,0,0),"mm")),
axis.title.y = element_text(size = base_size, angle = 90, vjust = 0.5, margin=unit(c(0,3,0,0),"mm")),
axis.ticks = element_line(size=base_line),
axis.ticks.length = unit(0.3, "lines"),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.border = element_rect(fill=NA, colour=NA, size=base_line),
panel.spacing = unit(1, "lines"),
legend.background = element_rect(fill="transparent", colour="transparent"),
legend.key = element_rect(fill="transparent", colour="transparent"),
legend.text= element_text(size=base_size),
strip.background = element_blank(),
strip.text.x = element_text(size = base_size * 0.8),
strip.text.y = element_text(size = base_size * 0.8, angle = -90),
strip.switch.pad.grid = unit(0, "mm"),
strip.switch.pad.wrap = unit(0, "mm"),
plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "cm"),
plot.title = element_text(size = base_size * 1.2),
plot.background = element_rect(colour = "transparent", fill="transparent", size=base_line)
)}Figure 1. Body parts of Megasoma gyas considered in the current study. From each body part, we obtained a measure of surface temperature over five minutes after animals were passively heated to 30 ºC and allowed to cool down to room temperature (~ 20 ºC). We used the change in surface temperature over time to obtain a measure of cooling rate for each body part.
Figure 2. Temporal changes in heat dissipation across body parts in Megasoma gyas.We used a water bath to heat beetles to 30 ºC and allowed them to dissipate heat to their surroundings (kept at 20 ºC) for 5 min. Each letter depicts a 1-min increment from the (A) start (t = 0 min) until the (F) completion (t = 5 min) of the experiment. The cephalic horn cooled faster than the other body parts considered in this study.
beetle<-read.csv(file = "cool-heat.csv", h = T, sep = ';')
beetle$invert<-beetle$cool.rate*-1 # obtaining positive cooling rate values
beetle2<-beetle[!(beetle$exp == "active.heating"),] # removing active heating experiments from data setbeetle.control<-subset(beetle2,exp=="control")
beetle.passive<-subset(beetle2,exp=="passive.heating")| ID | exp | body.part | cool.rate | mass | mass.cat | horn.ratio | invert | |
|---|---|---|---|---|---|---|---|---|
| 49 | C1 | control | cephal.horn | -1.0610 | 22.86 | high | 0 | 1.0610 |
| 50 | C1 | control | thorax.horn | -0.3569 | 22.86 | high | 0 | 0.3569 |
| 51 | C1 | control | elytra | 0.0000 | 22.86 | high | 0 | 0.0000 |
| 52 | C1 | control | abdomen | -0.3344 | 22.86 | high | 0 | 0.3344 |
| 53 | C2 | control | cephal.horn | -0.6196 | 22.86 | high | 0 | 0.6196 |
| 54 | C2 | control | thorax.horn | -0.7033 | 22.86 | high | 0 | 0.7033 |
| 55 | C2 | control | elytra | -0.3986 | 22.86 | high | 0 | 0.3986 |
| 56 | C2 | control | abdomen | -0.2050 | 22.86 | high | 0 | 0.2050 |
| 57 | C3 | control | cephal.horn | -0.5532 | 22.86 | high | 0 | 0.5532 |
| 58 | C3 | control | thorax.horn | -0.9675 | 22.86 | high | 0 | 0.9675 |
| 59 | C3 | control | elytra | -0.9305 | 22.86 | high | 0 | 0.9305 |
| 60 | C3 | control | abdomen | -1.0210 | 22.86 | high | 0 | 1.0210 |
| ID | exp | body.part | cool.rate | mass | mass.cat | horn.ratio | invert |
|---|---|---|---|---|---|---|---|
| 1 | passive.heating | cephal.horn | -0.3747 | 23.79 | high | 0.9540217 | 0.3747 |
| 1 | passive.heating | thorax.horn | -0.1000 | 23.79 | high | 0.9540217 | 0.1000 |
| 1 | passive.heating | elytra | -0.2329 | 23.79 | high | 0.9540217 | 0.2329 |
| 1 | passive.heating | abdomen | -0.0839 | 23.79 | high | 0.9540217 | 0.0839 |
| 2 | passive.heating | cephal.horn | -0.1852 | 22.90 | high | 0.6802248 | 0.1852 |
| 2 | passive.heating | thorax.horn | -0.0521 | 22.90 | high | 0.6802248 | 0.0521 |
| 2 | passive.heating | elytra | -0.0243 | 22.90 | high | 0.6802248 | 0.0243 |
| 2 | passive.heating | abdomen | -0.1726 | 22.90 | high | 0.6802248 | 0.1726 |
| 3 | passive.heating | cephal.horn | -0.3769 | 19.13 | low | 0.7269632 | 0.3769 |
| 3 | passive.heating | thorax.horn | -0.0227 | 19.13 | low | 0.7269632 | 0.0227 |
| 3 | passive.heating | elytra | -0.0412 | 19.13 | low | 0.7269632 | 0.0412 |
| 3 | passive.heating | abdomen | -0.1996 | 19.13 | low | 0.7269632 | 0.1996 |
| 4 | passive.heating | cephal.horn | -0.3623 | 23.32 | high | 0.6548562 | 0.3623 |
| 4 | passive.heating | thorax.horn | -0.0675 | 23.32 | high | 0.6548562 | 0.0675 |
| 4 | passive.heating | elytra | -0.0229 | 23.32 | high | 0.6548562 | 0.0229 |
| 4 | passive.heating | abdomen | -0.3104 | 23.32 | high | 0.6548562 | 0.3104 |
| 5 | passive.heating | cephal.horn | -1.2905 | 13.74 | low | 0.4223241 | 1.2905 |
| 5 | passive.heating | thorax.horn | -0.1514 | 13.74 | low | 0.4223241 | 0.1514 |
| 5 | passive.heating | elytra | -0.0918 | 13.74 | low | 0.4223241 | 0.0918 |
| 5 | passive.heating | abdomen | -0.1208 | 13.74 | low | 0.4223241 | 0.1208 |
| 6 | passive.heating | cephal.horn | -0.2828 | 22.86 | high | 0.7948937 | 0.2828 |
| 6 | passive.heating | thorax.horn | -0.0505 | 22.86 | high | 0.7948937 | 0.0505 |
| 6 | passive.heating | elytra | -0.0280 | 22.86 | high | 0.7948937 | 0.0280 |
| 6 | passive.heating | abdomen | -0.4571 | 22.86 | high | 0.7948937 | 0.4571 |
beetle.passive.mean<-beetle.passive %>%
group_by(body.part) %>%
summarise(avg = mean(invert), sd = sd(invert)); beetle.passive.mean## # A tibble: 4 × 3
## body.part avg sd
## <chr> <dbl> <dbl>
## 1 abdomen 0.224 0.138
## 2 cephal.horn 0.479 0.405
## 3 elytra 0.0735 0.0823
## 4 thorax.horn 0.0740 0.0455
beet.cont.mean<-beetle.control %>%
group_by(body.part) %>%
summarise(avg = mean(invert), sd = sd(invert)); beet.cont.mean## # A tibble: 4 × 3
## body.part avg sd
## <chr> <dbl> <dbl>
## 1 abdomen 0.520 0.439
## 2 cephal.horn 0.745 0.276
## 3 elytra 0.443 0.467
## 4 thorax.horn 0.676 0.306
We are accounting for body mass in the model
##
## Breusch-Godfrey test for serial correlation of order up to 8
##
## data: Residuals
## LM test = 7.0906, df = 8, p-value = 0.5269
##
## Call:
## lm(formula = invert ~ body.part + mass, data = beetle.passive)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26053 -0.08271 -0.02570 0.02655 0.65450
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.68076 0.26186 2.600 0.0176 *
## body.partcephal.horn 0.25467 0.11940 2.133 0.0462 *
## body.partelytra -0.15055 0.11940 -1.261 0.2226
## body.partthorax.horn -0.15003 0.11940 -1.257 0.2241
## mass -0.02179 0.01183 -1.842 0.0811 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2068 on 19 degrees of freedom
## Multiple R-squared: 0.4968, Adjusted R-squared: 0.3908
## F-statistic: 4.689 on 4 and 19 DF, p-value: 0.008396
## Df Sum Sq Mean Sq F value Pr(>F)
## body.part 3 0.6570 0.21901 5.121 0.00917 **
## mass 1 0.1452 0.14518 3.395 0.08107 .
## Residuals 19 0.8126 0.04277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
There is variation among body parts, with the cephalic horn cooling at a higher rate than other parts. Body mass did not mediate this effect.
| invert | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.68 | 0.13 – 1.23 | 0.018 |
| body part [cephal.horn] | 0.25 | 0.00 – 0.50 | 0.046 |
| body part [elytra] | -0.15 | -0.40 – 0.10 | 0.223 |
| body part [thorax.horn] | -0.15 | -0.40 – 0.10 | 0.224 |
| mass | -0.02 | -0.05 – 0.00 | 0.081 |
| Observations | 24 | ||
| R2 / R2 adjusted | 0.497 / 0.391 | ||
##
## Breusch-Godfrey test for serial correlation of order up to 8
##
## data: Residuals
## LM test = 4.9245, df = 8, p-value = 0.7656
##
## Call:
## lm(formula = invert ~ body.part + mass.cat, data = beetle.passive)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25640 -0.10392 -0.01390 0.01905 0.73749
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1869 0.0938 1.993 0.0608 .
## body.partcephal.horn 0.2547 0.1251 2.036 0.0559 .
## body.partelytra -0.1505 0.1251 -1.204 0.2435
## body.partthorax.horn -0.1500 0.1251 -1.200 0.2450
## mass.catlow 0.1114 0.0938 1.188 0.2496
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2166 on 19 degrees of freedom
## Multiple R-squared: 0.4479, Adjusted R-squared: 0.3316
## F-statistic: 3.853 on 4 and 19 DF, p-value: 0.01862
## Df Sum Sq Mean Sq F value Pr(>F)
## body.part 3 0.6570 0.21901 4.667 0.0132 *
## mass.cat 1 0.0662 0.06620 1.411 0.2496
## Residuals 19 0.8916 0.04693
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Fig1<-ggplot(beetle.passive, aes(x = mass, y = invert, colour=body.part, group=body.part, fill=body.part)) +
geom_segment(x=12,xend=24,y=beet.cont.mean$avg[1], yend=beet.cont.mean$avg[1],
linetype="dashed", color = "#EBC397", lwd=.5, alpha=.45)+
geom_segment(x=12,xend=24,y=beet.cont.mean$avg[2], yend=beet.cont.mean$avg[2],
linetype="dashed", color = "#A6E7FF", lwd=.5, alpha=.45)+
geom_segment(x=12,xend=24,y=beet.cont.mean$avg[3], yend=beet.cont.mean$avg[3],
linetype="dashed", color = "#2C809E", lwd=.5, alpha=.45)+
geom_segment(x=12,xend=24,y=beet.cont.mean$avg[4], yend=beet.cont.mean$avg[4],
linetype="dashed", color = "#9E5F1C", lwd=.5, alpha=.45)+
stat_smooth(formula = "y ~ x", method = "lm", fullrange = T, alpha=0.5,
show.legend = FALSE, se=F, lwd=.95) +
scale_colour_manual(values=c("#EBC397", "#A6E7FF","#2C809E", "#9E5F1C"))+
geom_point(shape = 21, colour = "black", size = 3, alpha = .9) +
scale_fill_manual(name = "", labels = c("Abdomen","Cephalic horn", "Thoracic horn","Elytra"),
values=c("#EBC397", "#A6E7FF","#2C809E", "#9E5F1C"))+
ylab("Cooling rate\n(slope)") +
xlab("Body mass (g)") +
theme(
# legend.position=c(0.8,0.85),
legend.title=element_text(size=11),
legend.text=element_text(size=10)) +
scale_y_continuous(limits=c(0,1.5), labels=seq(0,1.5,0.5))+
ggtheme(12,0.3)
# Beetle1 <- readPNG(paste0(figuresfolder, "/Beetle inset1.png"), native=T)
# Beetle2 <- readPNG(paste0(figuresfolder, "/Beetle inset2.png"), native=T)
Beetle3 <- readPNG(paste0(figuresfolder, "/Beetle inset3.png"), native=T)
Fig1Inset<- Fig1+
inset_element(p=Beetle3, 0.55, 0.7,1,.9); Fig1InsetFigure 3. Cooling rate as a function of body mass in Megasoma gyas. Body parts are color-coded, highlighting that heat dissipation occurs at different rates among body parts. The solid lines represent the expected relationship between cooling rates and body mass per body part. The dashed lines represent the cooling rates of different body parts measured from control specimens.
One individual (ID5) was much lighter than the others in our study. This individual had the cephalic horn with the highest cooling rate in our data set. To understand if the result we obtained was caused by a potential outlier, we will remove ID5 from our data set and re-run our analyses.
We are accounting for body mass in the model
##
## Breusch-Godfrey test for serial correlation of order up to 8
##
## data: Residuals
## LM test = 11.225, df = 8, p-value = 0.1893
##
## Call:
## lm(formula = invert ~ body.part + mass, data = beetle.passive.subset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16834 -0.04533 -0.01075 0.04340 0.20989
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.123462 0.298289 0.414 0.68480
## body.partcephal.horn 0.071660 0.062187 1.152 0.26722
## body.partelytra -0.174860 0.062187 -2.812 0.01314 *
## body.partthorax.horn -0.186160 0.062187 -2.994 0.00909 **
## mass 0.005413 0.013171 0.411 0.68688
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09833 on 15 degrees of freedom
## Multiple R-squared: 0.6318, Adjusted R-squared: 0.5336
## F-statistic: 6.434 on 4 and 15 DF, p-value: 0.003197
## Df Sum Sq Mean Sq F value Pr(>F)
## body.part 3 0.24717 0.08239 8.522 0.00153 **
## mass 1 0.00163 0.00163 0.169 0.68688
## Residuals 15 0.14502 0.00967
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Df Sum Sq Mean Sq F value Pr(>F)
## body.part 4 0.8415 0.21037 21.759 4.21e-06 ***
## mass 1 0.0016 0.00163 0.169 0.687
## Residuals 15 0.1450 0.00967
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Apparently, the cephalic horn and abdomen have cooling rates more similar to the control compared to the elytra and toracic horn. This makes sense, as flight muscles (which are known to play a role in insect thermoregulation) are concentrated in the region of the elytra/thorax.
##
## Breusch-Godfrey test for serial correlation of order up to 8
##
## data: Residuals
## LM test = 7.8528, df = 8, p-value = 0.448
##
## Call:
## lm(formula = invert ~ body.part + mass.cat, data = beetle.passive.subset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16389 -0.04585 -0.01375 0.04595 0.20931
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.24779 0.04546 5.450 6.7e-05 ***
## body.partcephal.horn 0.07166 0.06238 1.149 0.26861
## body.partelytra -0.17486 0.06238 -2.803 0.01337 *
## body.partthorax.horn -0.18616 0.06238 -2.985 0.00926 **
## mass.catlow -0.01535 0.05513 -0.278 0.78449
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09862 on 15 degrees of freedom
## Multiple R-squared: 0.6295, Adjusted R-squared: 0.5307
## F-statistic: 6.372 on 4 and 15 DF, p-value: 0.003335
## Df Sum Sq Mean Sq F value Pr(>F)
## body.part 3 0.24717 0.08239 8.471 0.00157 **
## mass.cat 1 0.00075 0.00075 0.078 0.78449
## Residuals 15 0.14590 0.00973
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Fig2<-ggplot(beetle.passive.subset, aes(x = mass, y = invert, colour=body.part, group=body.part, fill=body.part)) +
geom_segment(x=12,xend=24,y=beet.cont.mean$avg[1], yend=beet.cont.mean$avg[1],
linetype="dashed", color = "#EBC397", lwd=.5, alpha=.45)+
geom_segment(x=12,xend=24,y=beet.cont.mean$avg[2], yend=beet.cont.mean$avg[2],
linetype="dashed", color = "#A6E7FF", lwd=.5, alpha=.45)+
geom_segment(x=12,xend=24,y=beet.cont.mean$avg[3], yend=beet.cont.mean$avg[3],
linetype="dashed", color = "#2C809E", lwd=.5, alpha=.45)+
geom_segment(x=12,xend=24,y=beet.cont.mean$avg[4], yend=beet.cont.mean$avg[4],
linetype="dashed", color = "#9E5F1C", lwd=.5, alpha=.45)+
stat_smooth(formula = "y ~ x", method = "lm", fullrange = T, alpha=0.5,
show.legend = FALSE, se=F, lwd=.95) +
scale_colour_manual(values=c("#EBC397", "#A6E7FF","#2C809E", "#9E5F1C"))+
geom_point(shape = 21, colour = "black", size = 3, alpha = .9) +
scale_fill_manual(name = "", labels = c("Abdomen","Cephalic horn", "Thoracic horn","Elytra"),
values=c("#EBC397", "#A6E7FF","#2C809E", "#9E5F1C"))+
ylab("Cooling rate\n(slope)") +
xlab("Body mass (g)") +
theme(
# legend.position=c(0.8,0.85),
legend.title=element_text(size=11),
legend.text=element_text(size=10)) +
scale_y_continuous(limits=c(0,1.5), labels=seq(0,1.5,0.5))+
ggtheme(12,0.3);Fig2Notice the different x-axis dimensions.
R version 4.3.2 (2023-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8
attached base packages: stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: sjPlot(v.2.8.15), patchwork(v.1.2.0), png(v.0.1-8), pander(v.0.6.5), cowplot(v.1.1.2), kableExtra(v.1.3.4), forecast(v.8.21.1), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.2), readr(v.2.1.5), tidyr(v.1.3.0), tibble(v.3.2.1), ggplot2(v.3.5.1), tidyverse(v.2.0.0), lme4(v.1.1-35.1), Matrix(v.1.6-1.1), performance(v.0.10.8) and scales(v.1.3.0)
loaded via a namespace (and not attached): rlang(v.1.1.3), magrittr(v.2.0.3), tseries(v.0.10-55), compiler(v.4.3.2), mgcv(v.1.9-0), systemfonts(v.1.0.5), vctrs(v.0.6.5), quadprog(v.1.5-8), rvest(v.1.0.3), pkgconfig(v.2.0.3), fastmap(v.1.1.1), backports(v.1.4.1), labeling(v.0.4.3), effectsize(v.0.8.6), utf8(v.1.2.4), rmarkdown(v.2.25), tzdb(v.0.4.0), nloptr(v.2.0.3), xfun(v.0.41), cachem(v.1.0.8), jsonlite(v.1.8.8), highr(v.0.10), sjmisc(v.2.8.9), ggeffects(v.1.3.4), broom(v.1.0.5), parallel(v.4.3.2), R6(v.2.5.1), bslib(v.0.6.1), stringi(v.1.8.3), boot(v.1.3-28.1), lmtest(v.0.9-40), jquerylib(v.0.1.4), estimability(v.1.4.1), Rcpp(v.1.0.12), knitr(v.1.45), modelr(v.0.1.11), zoo(v.1.8-12), parameters(v.0.21.3), splines(v.4.3.2), nnet(v.7.3-19), timechange(v.0.2.0), tidyselect(v.1.2.0), rstudioapi(v.0.15.0), yaml(v.2.3.8), timeDate(v.4032.109), sjlabelled(v.1.2.0), curl(v.5.2.0), lattice(v.0.22-5), quantmod(v.0.4.25), withr(v.2.5.2), bayestestR(v.0.13.1), urca(v.1.3-3), coda(v.0.19-4.1), evaluate(v.0.23), xts(v.0.13.1), xml2(v.1.3.6), pillar(v.1.9.0), insight(v.0.19.7), generics(v.0.1.3), TTR(v.0.24.4), hms(v.1.1.3), munsell(v.0.5.0), minqa(v.1.2.6), xtable(v.1.8-4), glue(v.1.7.0), emmeans(v.1.9.0), tools(v.4.3.2), see(v.0.8.1), webshot(v.0.5.5), mvtnorm(v.1.2-4), grid(v.4.3.2), datawizard(v.0.9.1), colorspace(v.2.1-0), nlme(v.3.1-163), fracdiff(v.1.5-2), cli(v.3.6.2), fansi(v.1.0.6), viridisLite(v.0.4.2), svglite(v.2.1.3), sjstats(v.0.18.2), gtable(v.0.3.4), sass(v.0.4.8), digest(v.0.6.34), ggrepel(v.0.9.5), farver(v.2.1.1), htmltools(v.0.5.7), lifecycle(v.1.0.4), httr(v.1.4.7) and MASS(v.7.3-60)